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£H ■ Abstract 

a : 

The aim of this paper is to develop a crowd motion model designed to handle highly packed situations. The 
00 ■ model we propose rests on two principles: We first define a spontaneous velocity which corresponds to the 

velocity each individual would like to have in the absence of other people; The actual velocity is then computed 
as the projection of the spontaneous velocity onto the set of admissible velocities (i.e. velocities which do not 
violate the non-overlapping constraint). We describe here the underlying mathematical framework, and we 
■ explain how recent results by J.F. Edmond and L. Thibault on the sweeping process by uniformly prox-regular 

sets can be adapted to handle this situation in terms of well-posedness. We propose a numerical scheme for 
this contact dynamics model, based on a prediction-correction algorithm. Numerical illustrations are finally 
presented and discussed. 

Resume 
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!> ' Nous proposons un modele de mouvements de foule oriente vers la gestion de configurations tres denses. Ce 

modele repose sur deux principes: tout d'abord nous definissons une vitesse souhaitee correspondant a la 
vitesse que les individus aimeraient avoir en l'absence des autres; la vitesse reelle est alors obtenue comme 
projection de la vitesse souhaitee sur un ensemble de vitesses admissibles (i.e. qui respectent la contrainte de 
non-chevauchement). Nous decrivons le cadre mathematique sous-jacent et nous expliquons comment certains 
resultats de J.F. Edmond et L. Thibault sur les processus de rafle par des ensembles uniformement prox- 
reguliers peuvent etre utilises pour prouver le caractere bien pose de notre modele. Nous proposons un schema 
numerique pour ce modele de dynamique des contacts base sur un algorithme de type prediction-correction. 
Enfin des resultats numcriques sont presentes et commentes. 
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d 1 Introduction 



Walking behaviour of pedestrians has given rise to a large amount of empirical studies over the last decades. 
Qualitative data (preferences, walk tendencies) have been collected by Fruin [TS], Navin, Wheeler [37], Hender- 
son [20] and, more recently, by Weidmann [43] . From these observations, several strategies for crowd motion 
modelling have been proposed, and can be classified with respect to the way they handle people density (La- 
grangian description of individuals or macroscopic approach) , and to the nature of motion phenomena (deter- 
ministic or stochastic). Among discrete and stochastic models, let us mention Cellular Automata [21 [71 131>1 138j . 
models based on networks [TB] as route choice models [3J H] and queuing models [3TJ33]. In these models, each 
cell or node is either empty or occupied by a single person and people's motion always satisfies this rule. In 
cellular automata models, there are two manners of moving people during a time step. With the first one, po- 
sitions are updated one by one with a random order (Random Sequential Update). The second method consists 
of updating simultaneously all positions (Parallel Update). If several people want to reach the same cell, only 
one of them (randomly chosen) is allowed to move. In route choice models, people move on a network. Each 
model is based on a route choice set. Most choice set generation procedures are based on shortest route search 
and use shortest paths algorithms. Queuing models use Markov-chain models to describe how pedestrians move 
from one node of the network to another. 

In [19] , a microscopic model called social force model is presented. It describes crowd motion with a system 
of differential equations. The acceleration of an individual is obtained according to Newton's law. Several forces 
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are introduced as for example a term describing the acceleration towards the desired velocity or a repulsion 
force reflecting that a pedestrian tends to keep a certain distance from other people and obstacles. Moreover 
macroscopic models have been proposed. In |20j . pedestrian traffic dynamics is firstly compared with fluid 
dynamics. Some models [17l[20j[22] are based on gas-kinetic theory. Other models 23, 24, 25, 26 rest on a set 
of partial differential equations describing the conservation of flow equation. 

Several softwares have been developped: PedGo [21] . SimPed [11], Legion [40], Mipsim [22] or Exodus [16]. 
Some commonly observed collective patterns are now considered as standard benchmarks for those numerical 
simulations. Among these phenomena of self-organization, there is the formation of lanes formed naturally by 
people moving in opposite directions. In this way, strong interactions with oncoming pedestrians are reduced, 
and a higher walking speed is possible. Another phenomena is the formation of arches upstream the exit 
during the evacuation of a room. These patterns are recovered by CA-models [551 13S] and by the social force 
model QI1CEB]. 

The case of evacuation in emergency situations is of particular importance in terms of applications (obser- 
vance of security rules, computer-assisted design of public buildings, appropriate positionning of exit signs). 
Numerical simulations may allow to estimate evacuation time (to be compared for example with the duration 
of fire propagation) and also to predict areas where high density will appear. As pointed out by Helbing [18] . 
emergency situations do not fit into the standard framework of pedestrian traffic flow. When people stroll 
around without hurry, they tend to keep a certain distance from each other and from obstacles. In an emer- 
gency situation, the motion of individuals is governed by different rules. In particular, the contact with walls 
or other people is no longer avoided. Some strategies have been proposed to adapt social walk models to highly 
congested situation (see again [H]). We propose here an approach which relies on the very consideration that 
actual motion in emergency situations is governed by the opposition between achievement of individual satis- 
faction (people struggle to escape as quickly as possible, regardless of the global efficiency) and congestion. In 
particular, we aim at integrating the direct conflict between people in the model, in order to estimate in some 
way interaction forces between them, and therefore provide a way to estimate the local risk of casualties. 

The microscopic model we propose rests on two principles. On the one hand, each individual has a spon- 
taneous velocity that he would like to have in the absence of other people. On the other hand, the actual 
velocity must take into account congestion. Those two principles lead us to define the actual velocity field as 
the projection of the spontaneous velocity onto the set of admissible velocities (regarding the non-overlapping 
constraints). The flexibility of this model lies in its first point: every choice of spontaneous velocity can be 
made and so every existing model for predicting crowd motion can be integrated here. The key feature of the 
model is the second point which concerns handling of contacts. 

By specifying the link between these two velocities, the evolution problem takes the form of a first order 
differential inclusion. This type of evolution problem has been extensively studied in the 1970's, with the 
theory of maximal monotone operators (see e.g. [3]). A few years later, J.J. Moreau considered similar problems 
with time-dependent multivalued operator, namely sweeping processes by convex sets (see [35 ). Since then, 
important improvements have been developped by weakening the convexity assumption with the concept of 
prox-regularity. The well-poscdness of our evolution problem can be established by means of recent results of 
J.F. Edmond and L. Thibault [13j concerning sweeping processes by uniformly prox-regular sets. 

The paper is structured as follows: In Section 2, we present the model and establish its well-posedness; In 
Section 3, we propose a numerical scheme, and detail the overall solution method. Section 4 is devoted to some 
illustrations of the numerical algorithm. 

1 Modelling 

We consider N persons identified to rigid disks. For convenience, the disks are supposed here to have the same 
radius r. The centre of the i-th disk is denoted by q^ (see Fig. [1]). Since overlapping is forbidden, the vector of 
positions q = (qi,..,qjv) <= K 2Ar (equipped with the euclidean norm) is required to belong to the following set: 

Definition 1.1 (Set of feasible configurations) 

Q^qeK 2 ^, A,(q)>0 Vi<j}, 
where Dy(q) = |q^ — qj\ — 2r is the signed distance between disks i and j. 
We consider as given the vector of spontaneous velocities denoted by 

Utq^Ux^.^LMq^elR 2 ^ 
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Figure 1: Notations. 



Ui is the spontaneous velocity of individual i, which may depend on its own position (Uj = Uj(qj), see Section|4] 
for examples of such a situation), but also on other people's positions, that is why we keep here Uj = U,(q). 
To define the actual velocity, we introduce the following set: 

Definition 1.2 (Set of feasible velocities) 

C q = {veR 2W Vkj D«(q) = => Gy(q)-v>0}, 

with 

G„(q) = VAj(q) = (0,...,0,-e y (q),0,... I 0,e ii (q) ) 0,...,0) £ R 2N onde^q) = * ~ Qi 

The actual velocity field is defined as the feasible field which is the closest to U in the least square sense, which 
writes 

§ = P <- U <* (!) 
q(0) =q £Q, 

where Pc q denotes the euclidean projection onto the closed convex cone C q . 

Remark 1.3 Despite its formal simplicity, this model does not fit directly into a standard framework. Indeed 
the setCq does not continuously depend on q. If no contact holds, the velocity is not constrained andCq = M. 2N . 
With a single contact, the set C q becomes a half-space. 

2 Mathematical framework 
2.1 Reformulation 

Let us reformulate the problem by introducing A/q, the outward normal cone to the set of feasible configurations 
Q, which is defined as the polar cone of C q . 

Definition 2.1 (Outward normal cone) 

A/q = C° = {w £ R 2N , w • v < Vv e C q } • 



Remark 2.2 In Figure^ we represent the set Q C M. 2N which is defined as an intersection of convex sets' 
complements. In the case of a single contact (configuration qx), we remark that the cone Aqj is generated by 
the vector —G^qi) that is up to a constant, the outward normal vector to the domain D34 > 0. In the case of 
two or more contacts, the configuration q2 does not belong to a smooth surface and the cone Afq 2 (generated by 
— Gi2(q2) and — Gi3(q2)J generalizes somehow the notion of the outward normal direction. 

Thanks to Farkas' Lemma (see [8]), the outward normal cone can be expressed 

A/q = {- Ay G ii(<l) > ^3 > . Aj(q) > Ay = o} . (2) 

Let us recall the classical orthogonal decomposition of a Hilbert space as the sum of mutually polar cone 
(see [34]) : 

P Cq + PAA q = Id. (3) 
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Figure 2: Cones C q and A/q. 



Using this property, we get: 

^ = P Cq U(q) = U(q)-P A . q U(q). (4) 
Since Pjy U(q) E Aq, we obtain a new formulation for ([I]) 

§ + ^ 3U (* (5, 

q(0) = qo. 

The problem reads as a first order differential inclusion involving the multivalued operator Af. 

Remark 2.3 In the absence of contacts in the configuration q, the set of feasible velocities C q is equal to the 
whole space M. 2N , and consequently the outward normal cone TVq is reduced to {0}. In that case, the first relation 
of ® states that the actual velocity equals to the spontaneous velocity: 

If any contact exists, the differential inclusion means that the configuration q ; submitted to U(q), has to evolve 
while remaining in Q. 

Let us first study a special situation where standard theory can be applied. Consider N individuals in a corridor. 
In that case, as people cannot leap accross each other, it is natural to restrict the set of feasible configurations 
to one of its connected components: 

Q = {q = (<Jl, • • ■ ,qjv) 6 , q l+ i - q, > 2r}. 

In this very situation, as Q is closed and convex, the multivalued operator q i — ► A/" q identifies to the subdiffer- 
cntial of the indicatrix function of Q: 



dig (q) = {v, Iq (q) + (v, h) < I Q (q + h) Vh} , I Q (q) = 



if q e Q 

-oo if q ^ Q 



therefore q i — > Aq is maximal monotone. In that case, as soon as the spontaneous velocity is regular (say 
Lipschitz), standard theory (see e.g. Brezis [5]) ensures well-posedness. Yet, as illustrated in Figure^ Q is 
not convex in general and the operator q i — ► 7V q is not monotone. So we cannot apply the same arguments 
as in the case of a straight motion. By lack of convexity, the projection onto Q is not everywhere well-defined. 
However the set Q satisfies a weaker property in the sense that the projection onto Q is still well-defined in its 
neighbourhood. Indeed, Q is uniformly prox-regular, which is the suitable property to ensure well-posedness. 
Let us give some definitions to specify the general mathematical framework. 
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configuration q 



configuration q 
Figure 3: Lack of convexity. 



configuration q : 



q + q 




Figure 4: fy-prox-regular set. 



Definition 2.4 Let S be a closed subset of a Hilbert space H . 
We define the proximal normal cone to S at x by: 



wht 



N(S, x) = {veff, 3a > 0, x G P s (x + av)} , 



Ps(y) = {z e S, ds(y) = |y - z|}, ^ ds(y) = min|y - z| 



Following [10] . we define the concept of uniform prox- regularity as follows: 

Definition 2.5 Let S be a closed subset of a Hilbert space H. S is said rj-prox-regular if for all x £ dS and 
v G N(5, x), |v| — 1 we have: 

B(x + riv,n) H S = 0. 



In an euclidean space, S is 77-prox-regular if an external tangent ball with radius smaller than r\ can be rolled 
around it (see Fig 0|. Moreover, this definition ensures that the projection onto such a set is well-defined in its 
neighbourhood. The following remark will be useful later. 

Remark 2.6 If there exists a > satisfying x G Ps(x + cvv) then 

V/3>0, 13 < a, xGP S (x + /3v). 



Definition 2.7 The proximal subdifferential of function d$ at x is the set 

9 P ds(x) = |v G H, 3M, a>0, d s (y)-d s (x)+Af|y-x| 2 > (v,y-x>, VyGB(x,a)}. 

Let us specify the useful link between the previous subdifferential and the proximal normal cone, which is proved 
in 13121. 



Proposition 2.8 The following relation holds true: 

d p d s {x) = N P {S, x)n 5(0,1). 



Remark 2.9 A set C C H is convex if and only if it is oo-prox-regular. In this case N(C, x) = dlc(*-) for all 
xeC. 

We now come to the main result of this section. 

Theorem 2.10 Assume that U is Lipschitz and bounded. Then, for all T > and all qo S Q, the following 
problem 

ft +^U(q) 

q(o) = qo, 

has one and only one absolutely continuous solution q(-) over [0,T]. 

This well-posedness can be obtained by using results in [T3] [14] as soon as we prove that Q is uniformly prox- 
regular and that the set J\f„ identifies to the proximal normal cone to Q at q. This is the core of next subsection. 

Remark 2.11 It can be shown that the solution given by Theorem \ 2.10A satisfies the initial differential equa- 
tion (see 

2.2 Prox-regularity of Q 

Let us consider the set 

Q v ={qel M , A,(q) >0}. 

Proposition 2.12 Let S be a closed subset o/R™ whose boundary dS is an oriented C 2 hypersurface. For each 
x e dS, we denote by j/(x) the outward normal to S at x. Then, for each x € OS, the proximal normal cone to 
S at x is generated by ^(x) ; i.e. 

N(5,x) = R+z/(x). 



Proof: The proof is a straightforward computation (see [41]). □ 
We can also deduce the expression of the proximal normal cone to Qij . 

Corollary 2.13 For all q G Qij, 

N(Q <J -,q) = -R+G tf (q). 

By Definition 12. 5( the constant of prox-regularity equals to the largest radius of a "rolling external ball" . In 
order to estimate its radius, tools of differential geometry can be used. More precisely, to show that the set Qij 
is uniformly prox-regular, we can apply the following theorem, that is proved in [12] . 

Theorem 2.14 Let C be a closed convex subset o/R™ such that dC is an oriented C 2 hypersurface o/R™. We 
denote by z^c( x ) the outward normal to C at x and by pi(x), ..,p„_i(x) > the principal curvatures of C at x. 
We suppose that 

p= sup sup Pi(x) < oo. 

xedC l<i<n-l 

Then S = R" \ int(C) is a rj -prox-regular set with r\ = — . 

P 

Proposition 2.15 Qij is rj^-prox-regular with r]o — ry/2. 
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Proof: The set int(Qij) is obviously the complement of a convex set C which satisfies the assumptions of 
Theorem 12.141 The constant of prox-regularity of Qij can be obtained by calculating its principal curvatures, 
which are the eigenvalues of Weingarten endomorphism. Let q € dQij , the outward normal to C at q is equal 
to — f (q), where 

G tJ (q) (0,...,0,e I3 (q),0,...,0,-e I3 (q),0,...,0) 

Weingarten endomorphism is written as follows, for all tangent vectors h 6 Tq(dQij), 

W q (h) := -DKq)[h] = -— ^ r (o, . . • , 0, -P e x (h, - h, ; ), 0, . . . , 0, P c ^ (hj - h<), 0, . . . , o) , 

with 

P 6 _l Qxj - h^ = (hj - hi) - [(hj - hi) • ey]ey. 

After some computations, we deduce that the endomorphism W q has two eigenvalues, and y/2/\qj — q,-|, and 
the latter is equal to l/(rv / 2), which ends the proof. □ 
Now let us study the set of feasible configurations Q, that is the intersection of all sets Qij. We begin to 
determine its proximal normal cone. 

Proposition 2.16 For all q £ Q, N(Q,q) = ^N(Qy,q) = W q . 

Proof: The second equality follows from (|2|) and Proposition 12. 151 Let us prove the first one. If q S int(Q), 
then for each couple (i, j), q € int(Qy), which implies 

N(Q,q) = {0} = ^N(Q y ,q). 
We now consider q £ dQ and introduce the following set: 

Icontact = i < j, Ai(q) = 0} = {(i, j), i<j, q £ dQij}. (6) 

First, we check that N(Qy,q) C N(Q, q). Let belong to I con tact (otherwise the previous inclusion is 

obvious), we consider w £ N(<5y, q) \ {0} and we set v = w/|w|. By Proposition ^. 8[ v £ d p dq ij (q) and thus 

3M, a>0, d Qij (q) - d Qij (q) + M|q - q| 2 > v ■ (q - q), Vq £ B(q, a). 

Since dQ i3 -(q) = = dg(q) and rfQ y (q) < rfQ(q), it follows that 

3M, a>0, d Q (q)-d Q (q)+Af|q-q| 2 > v-(q-q), VqeB(q,a). 

Therefore v £ 9 p dg(q) and w £ N(Q, q). Consequently, for each couple £ Icontact, we obtain N(Qy , q) C N(Q, q) 
as required. We now want to prove 

It suffices to show that 

Vwi, w 2 £ N(Q,q)\{0}, w = wi +w 2 £ N(Q,q). 

Let Wi and w 2 belong to N(Q,q) \ {0}, we set w = Wj + w 2 , vi = wi/|wi| and v 2 = w 2 /|w 2 |. By Proposi- 
tion [2IS1 there exists Mi, M 2 > 0, ax, a 2 > such that 

d Q {c l )~d Q {q) + M 1 \q-q\ 2 > (vi,q-q>, VqeB(q, Ql ), 

dQ (q) - (q) + Af 2 1 q ~ q| 2 > (v 2 , q - q> , Vq G B (q, a 2 ) . 

So w = |wi|vi + |w 2 |v 2 and the vector v = w/|wi| + |w 2 | satisfies |v| < 1. Furthermore v = ivi + (1 — i)v 2 , 
where 

t = l w i| 

(|wi| + |w 2 |)' 

For a = min(ai, a 2 ) and M = tM\ + (1 — <)M 2 , the following relation holds 

d Q (q) - dq(q) + M|q - q| 2 > v • (q - q), Vq £ B(q, a). 
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Figure 5: Vanishing of the constant of prox- regularity. 



Hence v £ 3 p dg(q) and w £ N(Q, q). To conclude, it remains to check that 

N(Q,q) C^N^.q). 

By ([3]), any w £ N(Q, q) can be written w = v + z = P^w + Pc q w, with v_Lz. Suppose z ^ 0. Since 
w £ N(Q, q), there exists t > such that q G Pq{c[ + tw). Let 



s = min(t, e) with e = min 



(l,j)£lao n t a at \/2|z| 

by Remark 12.61 we know that q £ Pg(q + sw). Now set 



sw — sv 



and show that q £ Q. By convexity of Dij, we have 

Aj(q) > A,(q) + s Gy(q) • z, V(i,j). 
In addition, for (i, j) S I C ontact, it yields Gy(q) ■ z > 0, because z € C q . Consequently, 
V(i,i) £ Icontact, Aj(q) > Ai(q) + s Gij(q) • z = s G 8J (q) • z > 0. 

Furthermore, if (i, j) g Icontact, then s < ■ Hence 

V 2|z| 

Ai(q) > D«(q) + » Gii(q) ' z > Aj(q) - «V2|z| > 0. 

That is why q £ Q and dg(q+sw) < |q+sw— q| = s|v|. Yet |q+sw— q| = s|w| > s|v| because |w| 2 = |v| 2 + |z| 2 . 
Thus q ^ Pg(q + sw), which leads to a contradiction. In conclusion, z = and w = v £ 7V q = X3N(Qy,q), 
which completes the proof of the proposition. □ 
Now we want to show the uniform prox-regularity of Q. Since Q does not satisfy the same smoothness 
properties as Qij , the results of differential geometry cannot be applied. By Theorem 12. 141 if a set is the 
complement of a smooth convex set, then it is uniformly prox-regular. A natural question arises : Is the 
intersection of such sets (which is the case for Q) uniformly prox-regular with a constant depending only on the 
constants of prox-regularity of the smooth sets. From a general point of view, this is wrong as illustrated in 
Figure [5] Indeed, we have plotted in solid line the boundary of a set S which is the intersection of two identical 
disks' complements. This set is uniformly prox-regular but its constant of prox-regularity (equal to the radius 
of the disk plotted in dashed line) tends to zero when the disks' centres move away from each other. In this 
situation, the scalar product between normal vectors ni and n 2 (see Figure [6]) tends to -1. Thus, the constant 
of prox-regularity of S is also dependent on the angle between vectors ni and ri2. We now come to the main 
result of this subsection: the uniform prox-regularity of Q. This result rests on an inverse triangle inequality 
between vectors G„(q), which is based on angle estimates. Let us point out that we do not claim optimality 
of the constant r\ below. 

Proposition 2.17 Q is r] -prox-regular with 

n/2 1 



23AT ^SN 2 ' 
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Figure 6: Evolution of the angle between vectors ni and n 2 . 

Proof: We want to prove (cf. Proposition 12. 5p that there exists 77 > such that for all q G Q and for all 
veN(Q,q), 



v-(q-q) < ^|q-q| 2 , Vq G Q. 
By Proposition ^. 1 51 for all q G Qij and all w G N(Qy, q), we have 



(7) 



w • (q - q) < — |q - q| , Vq G Q lj . 
Inegality is obvious when v = 0. So we consider q G dQ and v G N(Q, q) \ {0}. By Proposition ^. 161 

(i,j)G.Icontact 



(8) 



We recall that Q C Qij so that by © we obtain 



(-S««G<i(q)) -(q-q)<S 



^^|q-q| 2 , VqeQ. 
2?7o 



The sum concerns only couples belonging to I CO ntact but for convenience, this point is omitted in the 

notation. As |Gj 3 -(q)| = y/2, we get 

v-(q-q) < ^(5>y) |q-q| 2 , v q e Q. 

To check Inequality (JT]), it suffices to find a constant ?/ > 0, independent from cnj and from q, satisfying 

1 1 



i.e. such that 



(E^)-y^<^|E^ G -(q) 



Finally, if we are able to exhibit 7 > verifying 



V2 



then Q will be 77-prox-regular with 



?7o r\/2 
»7 = — = 



7 7 

The problem takes the form of an inverse triangle inequality: 

E ai J I G « (q) I = E a w - 7 1 E a v Gi 3 (q) ■ 

The required result will follow as soon as we prove the main proposition stated below. 



□ 
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Proposition 2.18 (Inverse triangle inequality) 

There exists 7 > 1 such that for all q € Q, 



a y -|Gy(q)|<7 

contact 



contact 



whe 



Icontact — {(hj), i < j, -Dy (q) = 0} and are nonnegative reals. 
Constant 7 can be fixed as follows 



3N 



-1/2N 



Remark 2.19 TVoie the sign of coefficients onj. From a general point of view, this inequality is obviously wrong 
if these coefficients are just assumed real. Indeed, for N large enough, the cardinal of the set I con tact could 
be strictly larger than 2N , which induces a relation between vectors Gy(q) (see Fig. OH /or such a degenerate 
situation). 

The following elementary lemma asserts an inverse triangle inequality for two vectors. 

Lemma 2.20 Let u\ and U2 be two vectors ofM. 2N satisfying u\-u-i — cos0|ui||u2|, with cos9 > — 1. Then for 
all 



v>ve:= 



1 + cos t 



we have \u\\ + \u%\ < v\u\ +U2I. 



Proof of the inverse triangle inequality: We propose here a method based on angle estimates with vectors Gy (q) 
as pointed out in Figure [6] We use a recursive proof on the number of involved vectors. We are going to check 
that there exists S > 1 such that for all subset / C Icontact and for all otij > 0, 



Y a tJ \G tJ ( q )\<S^ 



(i,j')ei"c/co«toct 



{i,j)eICl C ontact 

Initialization: Suppose that the cardinality of I equals to 1, in other words, I = {(i, j)}. So we clearly have for 

ay|G«(q)| = KG y (q)| < *|a«G«(q)|. (9) 



all otij > and all S > 1 



Recursion assumption: 

If I J I — p, then we have for all otij > 



Y a y |G y (q)|<^ 

J <ZI contact 

Take a subset I C Icontact with |/| = p + 1. For any 



(i,j)£j<Zl contact 



(10) 



w = X <XijGij{q), 

with a.ij > 0, we choose (k, I) £ I and define J = I\ {(k, I)}, 

wi = 2J a u" G y (q) and w 2 = a fe iGfez(q). 



We need the following lemma which will be later proved. 
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-F k 

(a) Case 1 (b) Case 2a (c) Case 2b 



Lemma 2.21 J/wi 7^ 0, the following inequality holds 

ww / / 1 x 27VX 

> — k, urei/i k — 1 + — 



Wi Wj \ V 12 



Consequently, if wi =/= 0, from Lemma l2.20[ we deduce |wi | + |w2 1 < \ |wi + W2 1 (this inequality obviously 

\ 1 

holds for wi =0). By denoting S = J - > 1, we get 

|wi| + |w 2 | < <S|w|. (11) 
Applying recursion assumption (|10p and (|11[). we obtain 

Oij\Gij{q)\ < aw|Gw(q)| + ^ P |wi| < <5 P (|w 2 | + |w x |) < <5 p+1 |w|, 

which ends the proof of ([9]) by recursion. As \I CO ntact\ < 3iV, the inverse triangle inequality is checked with 
~f = 6 3N . □ 
Proof of Lemma T2. 2 II It suffices to deal with w 2 = G k i(q). By setting 



we have 



ctij if i < j 
ciji else, 



Wi = {F 1 ,F 2 ,...,F N ) where F p = ^/3 ip e lp . 

Thus, G K 2 can be interpreted as a pressure force exerted on the fc th person by its neighbours (different from 
the individual I). Similarly, — F k can be seen as a reaction force. We are looking for a lower bound of 

_ wi ■ w 2 _ -F k ■ e k i - F] ■ e ik 

|wi||w2 ' ~ ^VST^F' 

Case 1: -F k • e M > or -F t ■ e ik > 



Suppose that, for example (cf figure 7(a) I — Fk ■ eki > 0. Using \Fi ■ e/jt| < we get 

A > -Fi ■ e lk > -1 

In this case, k — 2 -1 / 2 . 

Case 2: —F k ■ &kl < and — i 7 } ■ e ik < 
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Case 2a: -F k ■ e kl > ~\F k \ or -Fj • e ;fe > ~\F t \ 

Suppose that, for example (cf Figure [7(b) [ ), — F k ■ e k i > — — \F k \. It can be shown that 

1 , — Ft • e k i , —Fi ■ e tk 

— - < — ■ and — ■ > — 1, 

which yields 

A . . > _ , - 

y/2 V 4 / 4\/2 

In this case k = 5/ (4v2). 



1 f 1 A 5 
A w > — I -- r-1 I = > -1. 



Case 2&: -F* • e fci < -^|F fe | and -Fj • e lk < -j|F;| (cf Figure [7(c) ) . 
We need the following lemma. 

Lemma 2.22 There exists k and I different from k and I verifying k ^ I and 

\F~ k \ > e\F k \, 

> e|F|, 

with e = 1/12 2N . 
We deduce that 



Therefore 



E l^l 2 ^ \ F «\ 2 + \ F '\ 2 + l^l 2 + \ F i\ 2 * (1 + e 2 ) [|F fc | 2 + |F;i| 2 ] 



ludes tl 

Proof of Lemma 12 221 We firstly consider 



In this case, k — , which concludes the proof of Lemma T2.21I □ 

Vl + e 2 



v k 

1=1 

where V k is the number of neighbours of individual k (individual I excepted) (V k < 5). As a consequence, 

v fe 

— F k ■ e k i — ^^^Pkj 0] i e kj _i ■ efci- 
i=i 

There exists ki e {jo,i,jo,2, -, ja,v k } (h ^ k,l) such that for all i e {1, ...,14} f3 kkl e kkl • e fci < P k j 0:i ekj ,i ' e fci- 
It is obvious that 

Pkki&kkx ■ Ski < ~Q^ k ' efci — — 24 I' 

In fact, individual fei is the neighbour who exerts the largest pressure force on person k. As illustrated in 
Figure individual k is between persons I and k\. 

If IFfcj | > — — | Ffc j , then we set k = k\. Else \F kl \ < -r^\F k \, and we produce the same reasoning with 
48 48 

—F kl — /3 klk e klk + 2^Pkiji,i e kih,ii 
i=i 

where V kl < 5. Thus, 

v kl 

—F kl ■ e k i = f3 klk e klk ■ e k i + '^^Pk 1 j 1:i ek 1 j 1 . i ■ 

i=l 



12 




Figure 7: Construction of sequence (fcj) 



Since -f3 klk e klk ■ e kl < ~^:\F k \ and -F kl ■ e u < \F kl \ < -r^\F k \, we obtain 



1 



As previously, there exists k% £ {ji,i,Ji,2j ■■■jJi,Vt 1 } (&2 ^ {^j^i}), such that 

1 



Pk 1 k 2 e k 1 k 2 • e fc2 < 



4 X 12 2 



1^1 



(Similarly, see Figure [3 individual ki is between persons ki and fc). 



1 / 1 



4 V 12 



|Ffe|, we set k = k%. Else, we continue by defining a sequence (ki) (cf Figure [7]) such that 
fco = k 



1 ( 1 



1 / 1 \ 1 



4 V 12 



0kik i+1 ekik i+1 ■ en < -j ( — J fj'"^'' 
It can be shown that fcj+i ^ {fco, fci, ..fci}. This construction ends at most in N — 2 steps: 

1 / 1 \ m 

3m < N - 1 satisfying |F fcm | > - f — j |F fc |. 



Finally we set 



Analoguously, we deal with Fj, by constructing a sequence verifying similar properties. We can check that 
^ (in proving that 

{fc , h, ..km} n {l 0l h, ..l p } = 0. 
The proof of Lemma \2. 221 is achieved by taking e = 1/12^. □ 



3 Numerical scheme 

3.1 Time-discretization scheme 

We present in this section a numerical scheme to approximate the solution to ((5J. The numerical scheme we 
propose is based on a first order expansion of the constraints expressed in terms of velocities. The time interval 
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is denoted by [0, T]. Let N 6 N*, h = T/N be the time step and t n — nh be the computational times. We 
denote by q™ the approximation of q(t n ). The next configuration is obtained as 

q n+l = q n + h 

where 

u" = P c h (U(q n )) with 

= {v G R 2N , Ai(q) + h Gytq) • v > Vz < j}. 
The scheme can be also interpreted in the following way. Let us introduce the set 

Q(q) ={qe M 2W , A,(q) + G« (q) • (q - q) > V* < j}, 

which can be seen as an inner convex approximation of Q with respect to q. Note that Q(q) is defined in such 
a way that Q is the union of all sets Q(q), q G Q- The scheme can be expressed in terms of position: 

q n+1 =PQ (q »)(q"+MJ(q n )). 

In this form it appears as a prediction-correction algorithm: predicted position vector q™ + /iU(q"), that may 
not be admissible, is projected onto the approximate set of feasible configurations. 

Remark 3.1 It is straightforward to check that 

n n+i _ n n 

h +N(Q(q"),q" +1 ) 3 U(q"), (12) 

so that the scheme can also be seen as a semi-implicit discretization of {5p, where N(Q(q"), q™ +1 ) approximates 
N(Q,q"). 

Convergence of this scheme shall be proven in a forthcoming paper. 
3.2 Numerical solutions 

In the model, the discrete actual velocity u™ is the projection of the spontaneous velocity onto the approximated 
set of feasible velocities. We propose here to solve this projection by a Uzawa algorithm (note that any algorithm 
could be used to perform this task). For convenience, explicit dependence of vectors and matrices upon the 
current configuration is omitted (e.g. U stands for U(q n ), for Dy(q™), etc. . . ). The actual velocity u solves 
the following minimization problem under constraints 

u = argmin |v — U| 2 . 
vec£ 

Uzawa algorithm is based on a reformulation of this minimization problem in a saddle-point form. We introduce 
the associated Lagrangian 

L (v, M ) - 1 1 v - U| 2 - V* ( D « + h G y ' v ) ' 

l<i<j<N 

and the following linear mapping 

N(N-l) 



B : R 2N - 

v h-> —h (Gij ■ v) 



l<3 



With these notations, the set Cq can be written 



i jv \ i\ — 1 ) 

= jv£M M ,V^ (R+)^^ , - ( D v + h G v ' v ) ^ 

l<i<j<W 

N(N — 1) 

v G R 2N , V/i G 2 , n ■ (By - D) < 
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Figure 8: A case of non-uniqueness for Kuhn- Tucker multipliers. 



where D — D(q) £ M. N ( N x '/ 2 is the vector of distances. The existence of a saddle-point 

(u,A) e M 2N x (R+)" 

for this problem is well-known (see e.g. [5]) and it is characterized by the next system: 

u + '.BA = U 
li ■ (Bu - D) < , V/i > 
A • (Bu-D) = 0. 

Uzawa algorithm produces two sequences (v fe ) G (R 2Ar ) N and (/i, fc ) e ^2 according to 

fi° = 

M fe+1 = n + ( M fc + P [Sv fc+1 -£>]), 

where 11+ is the euclidean projection onto the cone of vectors with nonnegative components (a simple cut-off in 
practice), and p > is a fixed parameter. The algorithm can be shown to converge as soon as < p < 2/\\B\\ 2 
(see |5J). More precisely, the sequence (v fc ) converges to u and it can be shown that the sequence (fi k ) tends 

to some A G (R + ) 5 such that (u, A) is a saddle-point of L. Notice that in general, the Kuhn- Tucker 
multiplier A is not unique as illustrated in Figure [8] In this case, the configuration of 14 people shows 29 
contacts, consequently matrix t B is not injcctive. 

Remark 3.2 (Link between local prox-regularity and speed of convergence for Uzawa algorithm) We denote by 
G the matrix whose columns are vectors Gij, where € Icontact (defined by fS)) J, and we introduce A = t GG. 
The size of this square matrix is equal to n conta ct which is the cardinal of Icontact- By inverse triangle inequality 
(see Proposition 1 2. 1 8\) . there exists a constant 7 such that for all A € (R+)"co»tod satisfying |A|i = 1, we have 



We define, for q € Q, a local parameter 7 q satisfying 



2 = t \ t GG\ = *AAA > 4- 

T 



2 

min l \A\ = — , 

IAh=l 7 2 



and ?7q = rv2/7q- Let us show that parameter ?7 q (setting a lower bound of the local prox-regularity of Q at 
point q) and the condition number of matrix A are closely related when A is non-singular. By denoting rj m i n 
the smallest eigenvalue of A, it follows that 

Vmin = min *AAA = min *AAA < min t \A\. 

A| 2 = l A| 2 >1 I*l2>i 

A>0 

Since for all A, |A|i < yfn contac t\ M2, we have 

min 'AAA < min 'AAA = n contact min 'AAA. 

|A| 2 >1 l*ll>V"contact l*ll>l 

A>0 A>0 A>0 
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Finally, 

Vmin < n con tact min 'AAA = n conta ct min l \A\ = 



l*U>i |Aix=i 72 

X>0 X>0 1 



6iV 

Furthermore, the condition number of matrix A equals to 

cond 2 (A) = UhWA- 1 ]] - Vmas 



2 



Since |Gy(q)| = \/2, we obtain \\A\\ 2 — rj max > 2, /i 



ence 



cond 2 (^) > > — > 



2 2 7 2 ^ 4 r 2 



r? mm " 67V " 6772 /V' 

which quantifies how the condition number of A varies with rjq. Since the matrix appearing in Uzawa algorithm 
is A — l GG, we expect that this algorithm converges less quickly for configurations with low local prox-regularity. 
In numerical simulations, we noticed indeed that solving the saddle-point problem requires more iterations in 
case of a jam. 



4 Numerical results 

In order to illustrate the contact model, we propose here an example of spontaneous velocity. The choice of 
the spontaneous velocity is important because this velocity reflects pedestrian behaviour. A lot of choices are 
obviously possible. The spontaneous velocity of an individual has to take into account obstacles in the room 
and specify how he wants to get around them. So this velocity depends on the room's geometry but it can 
be made dependent on other people positions too. Indeed, it is possible here to integrate individual strategies 
(deceleration or jam's avoiding). We refer the reader to [351 [331135] for other examples of spontaneous velocity. 
Here we restrict ourselves to simple behavourial model: people tend to optimize their own path, regardless of 
others. 

An example of spontaneous velocity 

We consider here the simplest choice for the spontaneous velocity. All the individuals have the same behaviour: 
they want to reach the exit by following the shortest path avoiding obstacles. Then, the spontaneous velocity's 
expression can be specified: 

U(q) = (Uo(qi), . . . , U (q w )) with U (x) - -S VD(x), 

where 2?(x) represents the geodesic distance between the position x and the nearest exit and s > denotes 
the speed. In order to compute T>, we have used the Fast Marching Method introduced by R. Kimmel and J. 
Sethian in [27] ■ In this method, the value of D is computed at each point of a grid. The value at the exit's 
nodes is set to zero. Then, the values of the distance at the other points is computed step by step so that a 
discrete version of |V2?| = 1 is satisfied. Moreover, the distance at the nodes situated in the obstacles is fixed 
to a large value, which prevents the shortest path from going across them. In Figure we have considered a 
room with 5 obstacles and the exit is situated to the left. We note that by following the built velocity field, 
people are going to avoid obstacles. 

Our aim is to simulate evacuation of any building consisting of several floors. We have chosen an object 
oriented programming method and we have implemented this Fast Marching Method in a C++ code. Let us 
detail this code. On each floor, the spontaneous velocity is directed by the shortest path avoiding obstacles to 
the nearest exit or stairwell. In the stairs, people just want to go down. We have integrated this spontaneous 
velocity in the C++ code SCoPI: Simulations of Collections of Interacting Particles developped by A. Lefebvre 
(see [2j2[3n]). This code allows us to compute the actual velocity as the projection of the spontaneous velocity 
as described in Section [3J 



l(i 




Remark 4.1 Notice that the velocity field produced by this strategy is not continuous as soon as the room is 
not convex, which rules out Theorem \2.1(A This lack of regularity is not important in practical applications : 
the places at which it occurs (in particular upstream obstacles) are emptied after a few moments. The main 
consequence is the discontinuity of the future configurations with respect to initial data, which is not surprising 
from a modelling standpoint. 

We propose to illustrate the behaviour of the algorithm in two situations. The first one corresponds to a 
many- individual evacuation from a square room through a single exit, the second one illustrate the capability of 
the approach to handle complicated geometries. For these two experiments, it will be noticed that the contacts 
between the individuals and the obstacles have to be handled (as the contacts between people). Even if an 
individual want to avoid an obstacle, he can be pushed on it by people behind them. 

Simple evacuation 

We consider the situation of 1000 people which are randomly distributed over a square room. The spontaneous 
velocity field corresponds straight pathlines towards the exit at constant speed. As the field has a negative 
divergence, it tends to increase the local density, so that congestion is rapidly reached in the neighbourhood 
of the exit, and the congestion front propagated upstream as long as it is feeded by incoming people. In 
Figure [TT] we represented the current configuration and the corresponding network of interaction pressures: for 
any couple of disks in contact, we represent the segment between centers, having its color (from white to black) 
depend upon the (positive) Kuhn- Tucker multiplier which handles the corresponding contraint. We recover the 
apparition of arches upstream the exit. The Kuhn- Tucker multipliers Xij quantify the way U, the spontaneous 
velocity field, does not fit the constraints, and as such they can be interpreted in terms of pressures undergone 
by individuals. Although it would be presumptuous at this stage to assimilate Xij to an actual measure of the 
discomfort experienced by persons i and j, it is obvious that high values for those Kuhn- Tucker multipliers can 
be expected on zones where people are likely to be crushed. 

Complex geometry 

In the second example we consider the evacuation of a floor through exit stairs. A zoom on the geometry 
near the exit (together with the isovalues of the geodesic distance function, on which the spontaneous velocity 
is built) is represented on Figure flOl Figure fl2l corresponds to snapshots at times 0s, 5s, lis, 16s, 41s and 
75s. Disks are colored according to their initial geodesic distance to the exit. Note that initial ordering is not 
preserved during the evacuation. Notice also how a jam forms between snapshots 2 and 3 in the room located 
on the left hand side. This jam decreases significantly the rate at which people exit the room, but it disappears 
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Figure 10: Geometry and isovalues for the geodesic distance. 



eventually. The final evacuation time is 109s, to be compared to 48s which corresponds to the evacuation time 
without congestion. 
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